#if !defined( MATRIX_H)
#define  MATRIX_H
#include <algorithm>
#include <memory>
#include <iostream>
#include <vector>
#include <list>
#include <string>
#include <cassert>
#include <functional>
#include "stlhelper.h"
#include <bits/stl_function.h>
#include <boost/scoped_array.hpp>

#include "matrix_traits.h"
#include "matrix_common.h"
namespace bcpp_matrix
{


template <typename T> class t_dense_matrix
{
private:
	boost::scoped_array<T>					raw_data;
	std::vector<T*>							data;
	void copy_from_same_size(const t_dense_matrix& other)
	{
		// copy from other to me
		for (unsigned j = 0; j < data.size(); ++j)
			std::copy(	other.data[j].get(), other.data[j].get() + data.size(),
						data[j].get());
	}
	T*										raw_begin() {return raw_data.get();}
	T*										raw_end()  {return raw_data.get() + size() * size();}
	const T*								raw_begin() const  {return raw_data.get();}
	const T*								raw_end() const   {return raw_data.get() + size() * size();}
public:
	unsigned size() const {return data.size();}
	void fill(T t)
	{
		std::fill(raw_begin(), raw_end(), t);
	}
	bool operator != (t_dense_matrix& other) const
	{
		return ! operator == (other);
	}
	bool operator == (t_dense_matrix& other) const
	{
		if (size() != other.size())
			return false;
		return std::equal(raw_begin(),  raw_end(),  other.raw_begin());
	}

	template <typename OP>
	t_dense_matrix& apply_operator(const t_dense_matrix& other, OP op)
	{
		if (size() == other.size())
			std::transform(raw_begin(), raw_end(), other.raw_begin(),  raw_begin(), op);
		else
		{
			unsigned sz = std::min(size(),  other.size());
			for (unsigned i = 0; i < sz; ++i)
				std::transform(data[i], data[i] + sz, other.data[i],  data[i], op);

		}
		return *this;
	}

	t_dense_matrix& operator *= (const t_dense_matrix& other)
	{
		return apply_operator(other, std::multiplies<T>());
	}
	t_dense_matrix& operator /= (const t_dense_matrix& other)
	{
		return apply_operator(other, std::divides<T>());
	}
	t_dense_matrix& operator -= (const t_dense_matrix& other)
	{
		return apply_operator(other, std::minus<T>());
	}
	t_dense_matrix& operator += (const t_dense_matrix& other)
	{
		return apply_operator(other, std::plus<T>());
	}
	t_dense_matrix& operator |= (const t_dense_matrix& other)
	{
		return apply_operator(other, bit_or<T>());
	}
	t_dense_matrix& operator &= (const t_dense_matrix& other)
	{
		return apply_operator(other, bit_and<T>());
	}


	void swap(t_dense_matrix& other)
	{
		raw_data.swap(other.raw_data);
		data.swap(other.data);
	}

	//
	// constructor
	//
    typedef T element_type;
	t_dense_matrix(){}
	t_dense_matrix(unsigned sz)
	{
		allocate(sz);
	}
	t_dense_matrix(unsigned sz, T value)
	{
		allocate(sz, value);
	}
	t_dense_matrix(const t_dense_matrix& other)
	{
		allocate(other.size());
		copy_from(other);
	}
	t_dense_matrix& operator = (const t_dense_matrix& other)
	{
		if (&other == this)
			return *this;

		// reallocate me if necessary
		if (size() != other.size())
			allocate(other.size());
		copy_from(other);

		return *this;
	}
	void min(const t_dense_matrix& other)
	{
		if (size() == other.size())
		{
			// same size: use raw data
			for (int i = size() * size(); i >= 0; --i)
				raw_data[i] = std::min(raw_data[i], other.raw_data[i]);
		}
		else
		{
			unsigned sz = std::min(size(), other.size());
			for (unsigned i = 0; i < sz; ++i)
				for (unsigned j = 0; j < sz; ++j)
					data[i][j] = std::min(data[i][j], other.data[i][j]);
		}
	}

	void max(const t_dense_matrix& other)
	{
		if (size() == other.size())
		{
			// same size: use raw data
			for (int i = size() * size(); i >= 0; --i)
				raw_data[i] = std::max(raw_data[i], other.raw_data[i]);
		}
		else
		{
			unsigned sz = std::min(size(), other.size());
			for (unsigned i = 0; i < sz; ++i)
				for (unsigned j = 0; j < sz; ++j)
					data[i][j] = std::max(data[i][j], other.data[i][j]);
		}
	}

	void copy_from(const t_dense_matrix& other)
	{
		partial_copy_from(other, other.size());
	}
	void partial_copy_from(const t_dense_matrix& other, unsigned sz)
	{
		assert(sz <= other.size());
		assert(sz <= size());
		if (sz == size())
		{
			std::copy (other.raw_begin(),  other.raw_end(), raw_begin());
		}
		else
		{
			for (unsigned i = 0; i < sz; i++)
				std::copy (other.data[i], other.data[i] + sz, data[i]);
		}

	}

	T& 			operator()(unsigned i,  unsigned j) { return data[i][j];	}
	const T&	operator()(unsigned i,  unsigned j) const { return data[i][j];	}
	T* operator[](unsigned i)  {assert(i < size()); return data[i];}
	const T* operator[](unsigned i) const {assert(i < size()); return data[i];}
	void allocate(unsigned sz, T value)
	{
		allocate(sz);

		// fill with default value
		fill(value);
	}
	void allocate(unsigned sz)
	{
		// reallocate me if necessary
        //std::cerr << sz * sz * sizeof(T) << "\n";
		if (size() != sz)
		{
			//std::cerr << "allocate\n";
			raw_data.reset(new T [sz * sz]);
			data.clear();
			for (T* i = raw_begin(), *end = raw_begin() + sz * sz; i != end; i += sz)
				data.push_back(i);
		}
	}

	t_dense_matrix&	sub_matrix(
									const std::vector<unsigned>&	indices,
									t_dense_matrix&				new_matrix) const
	{
		if (!indices.size())
			return new_matrix;
		typedef std::vector<unsigned>::const_iterator t_uiter;
		new_matrix.allocate(indices.size());

		for (unsigned i = 0; i < indices.size(); ++i)
		{
			for (unsigned j = 0; j < indices.size(); ++j)
			{
				assert(indices[i] < size());
				assert(indices[j] < size());
				new_matrix[i][j] = data[indices[i]][indices[j]];
			}
		}
		return new_matrix;
	}


};
template <typename T> inline void
swap(t_dense_matrix<T>& lhs, t_dense_matrix<T>& rhs) {lhs.swap(rhs);}

typedef t_dense_matrix<float> t_dense_matrix_f;
typedef t_dense_matrix<char> t_dense_matrix_bool;

template <typename T> inline
const t_dense_matrix<T> operator - (const t_dense_matrix<T>& lhs, const t_dense_matrix<T>& rhs)
{return t_dense_matrix<T>(lhs) -= rhs;}
template <typename T> inline
const t_dense_matrix<T> operator + (const t_dense_matrix<T>& lhs, const t_dense_matrix<T>& rhs)
{return t_dense_matrix<T>(lhs) += rhs;}
template <typename T> inline
const t_dense_matrix<T> operator | (const t_dense_matrix<T>& lhs, const t_dense_matrix<T>& rhs)
{return t_dense_matrix<T>(lhs) |= rhs;}
template <typename T> inline
const t_dense_matrix<T> operator & (const t_dense_matrix<T>& lhs, const t_dense_matrix<T>& rhs)
{return t_dense_matrix<T>(lhs) &= rhs;}
template <typename T> inline
const t_dense_matrix<T> operator / (const t_dense_matrix<T>& lhs, const t_dense_matrix<T>& rhs)
{return t_dense_matrix<T>(lhs) /= rhs;}
template <typename T> inline
const t_dense_matrix<T> operator * (const t_dense_matrix<T>& lhs, const t_dense_matrix<T>& rhs)
{return t_dense_matrix<T>(lhs) *= rhs;}



template <typename T, typename STRM> inline
void dump_to_stream(STRM& os, const t_dense_matrix<T>& matrix)
{
//	matrix in square
	for (unsigned i = 0, sz = matrix.size(); i < sz; ++i)
	{
		for (unsigned j = 0; j < sz; ++j)
			os << "\t" << matrix[i][j];
		os << "\n";
	}
}

template <typename T, typename STRM> inline
void dump_to_stream( STRM& os, const t_dense_matrix<T>& matrix, const std::vector<std::string>& names)
{
//	matrix in square
	for (unsigned i = 0, sz = matrix.size(); i < sz; ++i)
		os << "\t" << names[i];
	os << "\n";
	for (unsigned i = 0, sz = matrix.size(); i < sz; ++i)
	{
		os << names[i];
		for (unsigned j = 0; j < sz; ++j)
			os << "\t" << matrix[i][j];
		os << "\n";
	}
}


//________________________________________________________________________________________


// Holding column names, matrix names, nulls


//________________________________________________________________________________________
template <typename T> 
struct t_matrix
{
private:
	unsigned					cnt_nulls;
public:
	t_dense_matrix<T>			matrix;
	std::string					matrix_name;
	std::vector<std::string>	names;
	T	 						null_value;
	bool						is_missing(unsigned i, unsigned j) const
								{
									return null_value == matrix[i][j];
								}
	void						erase_element(unsigned i,  unsigned j)
										{matrix[i][j] = null_value;}
	unsigned					size() const {return matrix.size();}
	t_matrix<T>(T null_value_ = t_matrix_traits<T>::binary_null()):
											cnt_nulls(0),
											null_value(null_value_){}
	~t_matrix(){}
	// does not assume symmetrical matrix
	unsigned					get_count_filled() const
										{ return get_count_pairs() - get_count_missing();}
	unsigned					get_count_pairs() const
										{ return size () * size();}
	unsigned					get_count_missing() const
										{ return cnt_nulls;}
	void allocate(unsigned size, T null_value_);
	void allocate(unsigned size);
	void get_single_linkage_clusters(std::vector<t_matrix<T> >& data) const;

	void count_missing();
	bool operator == (t_matrix<T>& other) const
	{
		return  matrix_name == other.matrix_name &&
				matrix == other.matrix &&
				names == other.names &&
				null_value == other.null_value;
	}
	void swap(t_dense_matrix<T>& other, const std::vector<std::string>& names_, 
				const T& null_  = t_matrix_traits<T>::binary_null())
	{
		names = names_;
		matrix.swap(other);
		cnt_nulls = 0;
		null_value = null_;
	}


	void swap(t_matrix<T>& other)
	{
		matrix.swap(other.matrix);
		names.swap(other.names);
		std::swap(null_value, other.null_value);
		std::swap(cnt_nulls, other.cnt_nulls);
	}

	void reduce_to_subset(const std::vector<unsigned>& order_indices);
	void subset(const std::vector<unsigned>& indices, t_matrix<T>& new_data) const;
	void sort_by_name_alpha();
};
template <typename T> 
inline void swap(t_matrix<T>& lhs, t_matrix<T> rhs) {lhs.swap(rhs);}


typedef t_matrix<float> t_matrix_f;


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	single linkage clusters

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

template <class T>
void t_matrix<T>::get_single_linkage_clusters(std::vector<t_matrix<T> >& data) const
{

	//
	//	separate matrix using single linkage clustering
	//
	std::vector<std::vector<unsigned> > clusters;
	private_details::get_single_linkage_clusters_indices(*this, clusters);
	if (clusters.size() == 1)
	{
		data.push_back(t_matrix<T>(null_value));
		data.back() = *this;
		return;
	}

	//
	//	for each cluster
	//
	for (unsigned i = 0; i < clusters.size(); ++i)
	{
		data.push_back(t_matrix<T>(null_value));

		// divide up data
		matrix.sub_matrix(clusters[i], data.back().matrix);
		data.back().count_missing();
		// divide up names
		private_details::extract_clustered_names(names , clusters[i], data.back().names);
		// data' name is my name
		data.back().matrix_name = matrix_name;
	}
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	t_matrix_io functions

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template <class T>
void t_matrix<T>::allocate(unsigned size, T binary_null)
{
	null_value = binary_null;
	allocate(size);
}

template <class T>
void t_matrix<T>::allocate(unsigned size)
{
	matrix.allocate(size, null_value);
	names.assign(size, "");
}

template <typename T> 
void t_matrix<T>::count_missing()
{
	cnt_nulls = 0;
	for (unsigned i = 1, sz = size(); i < sz; ++i)
		for (unsigned j = 0; j < i; ++j)
			if (matrix[i][j] == null_value)
				++cnt_nulls;
};

template <typename T> 
void t_matrix<T>::reduce_to_subset(const std::vector<unsigned>& order_indices)
{
	t_matrix<T> new_data;
	subset(order_indices, new_data);
	swap(new_data);
}

template <typename T> 
void t_matrix<T>::subset(const std::vector<unsigned>& indices, t_matrix<T>& new_data) const
{
	// copy values to new data
	new_data.null_value	= null_value;
	new_data.matrix_name= matrix_name;

	// get data subset
	matrix.sub_matrix(indices, new_data.matrix);


	// get names subset
	for (unsigned i = 0; i < indices.size(); ++i)
		new_data.names.push_back(names[indices[i]]);

	// update missing count
	new_data.count_missing();
}





//} // private_details

template <typename T> 
void t_matrix<T>::sort_by_name_alpha()
{
	std::vector<unsigned> indices;
	sort_indices(names, indices);
	reduce_to_subset(indices);
}


} // bcpp_matrix



#endif	//MATRIX_H

